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This article is based on a talk given by one of us (EVI) at the conference "StatPhys- Taipei- 1997" . 
It overviews the exact results in the theory of the sandpile model and discusses shortly yet unsolved 
problem of calculation of avalanche distribution exponents. The key ingredients include the analogy 
with the critical reaction-diffusion system, the spanning tree representation of height configurations 
and the decomposition of the avalanche process into waves of topplings. 



PACS number(s): 05.40. +j 
I. INTRODUCTION 

The concept of Self-Organized Criticality one of 
the most fundamental concepts of modern physics of non- 
equilibrium phenomena; actually has a long history, dat- 
ing back about a half of the century, when Kolmogorov 
put forward his seminal theory of isotropic and homoge- 
neous turbulence 0. The cornerstone of his very simple, 
but surprisingly robust dimensional analysis, was the as- 
sumption that the fluid being driven by a random force 
evolves to a stationary state where the velocity correla- 
tion functions are universal and obey power laws, pro- 
vided that the viscosity of the fluid is small enough to 
ensure the existence of the so-called inertial range of 
scales. It was assumed that within this inertial interval 
the energy is being transferring from large eddies to small 
ones without any dissipation whatsoever. Such a non- 
equilibrium but stationary state, with a constant flux of 
conserved quantities, is usually called the Flux State. 

Later on, Kardar, Parisi and Zhang j^] retraced very 
nearly the same path on an absolutely different ground. 
Namely, they considered the simple type of growth exem- 
plified by a vapor-deposition process, where the growth 
rate is locally determined by the flux of particles arriv- 
ing ballistically at the surface. They found that this sys- 
tem evolves to the non-equilibrium stationary state where 
spatial and temporal fluctuations of the growing surface 
obey power laws. 

After a careful analysis of all these rather phenomeno- 
logical theories, Bak, Tang and Wiesenfeld (|] have 
squeezed the concept of Self-Organized Criticality out in 
its most purified form. They introduced a cellular au- 
tomaton now commonly known as " sandpile" because of 
the crude analogy between its dynamical rules and the 
way sand topples when building a real sand pile. 

The formulation of this model is given in terms of in- 
teger height variables Zi at each site of the square lattice 
C. Particles are added randomly and the addition of a 
particle increases the height at that site by one. If this 
height exceeds the critical value z c = 4, then the site 
topples, on toppling its height decreases by 4 and the 
heights at each of its nearest neighbors increases by 1. 



These may become unstable in their turn and the relax- 
ation process continues. This chain reaction propagates 
up to the moment when all sites become stable again. 
One assumes the updating to be done concurrently, with 
all sites updated simultaneously. Open boundary condi- 
tions are usually assumed, so that the toppled boundary 
site gives one particle to each of its three neighbors while 
one grain drops out of the system. 

This system also evolves stochastically into a critical 
state with a constant flux of particles on which it exhibits 
properties similar to that of a second order phase tran- 
sition Q. It lacks therein any characteristic length or 
time scale and obeys power-law distributions. The crit- 
ical properties of the flux state are independent of the 
initial configuration of the system and are determined 
completely by the flux rate. Unlike ordinary critical phe- 
nomena, no fine tuning of any other control parameters 
is necessary to arrive at this state. 

II. MEAN-FIELD APPROXIMATION 

A. Flux State 

For our derivation of the mean-field equations it will be 
even more convenient to generalize the dynamical rules of 
the model [|| . Namely, let us suppose that at each site of 
the lattice one of the species A, B, C or D is living. These 
species represent respectively four stable heights of origi- 
nal model. Then, due to the external driving force, some 
new particles (we will call them if) are added randomly 
into the system initiating avalanches. Actually, if are 
the only mobile species in this model and the avalanche 
process can be described in terms of their propagation 
through the correlated media of immobile species A, B, 
C and D. In this propagation the t^-species mutate all 
others in the sites they have visited and topple the criti- 
cal ones according to the following rules 

A + ^^B, 
B + <p->C, 

C + ¥>->D, (1) 
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Here ip and if denote the particles obtained by the site 
and the particles transferred to the neighboring sites after 
toppling, respectively. These processes can formally be 
reinterpreted as an irreversible chemical reaction which 
takes place at each site of the lattice or, in other words, 
as a branching process which is characterized by the so- 
called branching probabilities 

p = (pi,p 2 ,P3,Pi), Pi +P2 +Pa + P& = 1- (2) 

This generalized model also describes the critical flux 
state where the matter is conserved. It (obviously) corre- 
sponds to the total conservation of particles in the origi- 
nal formulation of the sandpile model. The latter can be 
reproduced exactly if we choose p = (0, 0, 0, 1). 

The simplest possible description of the critical flux 
state can be achieved within the so-called mean-field ap- 
proximation. There are a lot of different methods to cal- 
culate mean- field critical exponents All of them 
are simple enough and have many features in common 
that, when applied to our model of chemical reactions, 
can be summarized as follows. 

At first, we introduce the concentrations of species A, 
B, C and D 

n= (iiA,n B ,Ji c ,n D ), n A + n B + n c + riD = 1- (3) 

The normalization condition comes from the constraint 
that we always have only one species in each given site. 

Then, following the standard prescriptions of chemical 
physics, we can write kinetic equations corresponding to 
this scheme of chemical reactions 
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where p — p\ + 2p 2 + 3p3 + 4j»4 is equal to the average 
number of particles ip leaving the cell on toppling and r 
is the position vector of the site in the 2D space. The 
noise term ry(r,i), being non-negative, mimics the ran- 
dom addition of particles to the system. The diffusion 
term V 2 (n (/J riD) describes the transfer of particles into 
the neighboring cells, and the diffusion coefficient v for 
the discrete Laplacian on the square lattice is equal to 

The physical meaning of these equations is transpar- 
ent. When the concentration n v of species if is equal to 
zero, all toppling processes die. Then, due to the noise 
term ri{v,t), the particles are added randomly into the 



system initiating a branching process directed towards 
the open boundary of the system. This process mutates 
species in the cells it has visited and topples the criti- 
cal ones. Finally, the system will reach the steady state 
where the probability that the activity will die is on av- 
erage balanced by the probability that the activity will 
branch. Thus, the chain reaction maintains this station- 
ary state and all further avalanches cannot change the 
concentrations of species A, B, C, and D. Therefore, the 
steady state is characterized by the conditions that 



h B = h c = ho = 



(5) 



and Eqs. ([la-4d) lead to the following relations between 
concentrations of species n at the stationary state and 
branching probabilities p 



n l = Pi I Pi 

n B = (P3+Pi)/P, 

n C = (P2+P3+Pi)/P, 

n B = (Pi +P2+P3+ Pi)/P = l/P 



(6a) 
(6b) 
(6c) 
(6d) 



Here, within the mean-field approximation, we obviously 
neglected all spatial and temporal fluctuations of concen- 
trations n. 



B. Branching Process 

The next question is how the avalanche process can be 
described within the same approximation. It can natu- 
rally be represented as the critical branching process. To 
describe it we first subdivide the total number of sites 
toppled in the avalanche into the four terms: N A , Nb, Nc 
and TVd . These correspond respectively to the numbers of 
sites that contain species A, B, C and D after toppling in 
the avalanche. Their sum N = Na + N b + Nc + N-£> corre- 
sponds to the total number of topplings in the avalanche 
at the moment t. 

Then, we introduce the probability 



P t (N A ,NB,N c ,N u ), 



(7) 



of having such numbers after t time steps of avalanche 
process and the corresponding generating function 

Gt(^A, Ab, Ac, Ad) = 

Pt (N A , N B , A c , N u ) A^ A A^ B A^ c A£ d . (8) 

N 

Initially, at the moment t = 0, the only nonzero ele- 
ment of the probability distribution is Po(0, 0,0,0) = 1, 
which correspond to the value Go = 1 of the generating 
function. 

On the next time step, after the toppling of the first 
site, all nonzero probabilities are follows: 
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Pl(l, 0,0,0) = P 4UB, 

P 1 (0,l,0,0)=p 3 n D + n A , 
P 1 (0,0,l,0)=p 2 n D + n B , 
Pi (0,0, 0,1) = Pl n B + n c , 

with corresponding generating function 

Gi(Aa, Ab, Ac, Ad) = 

(n D P4)AA + (noP3 + n A )X B + 
(nr>P2 + n B )Xc + ("dPi + ™c)Ad, 



(9) 
(10) 

(11) 
(12) 



(13) 
(14) 
(15) 



Now, it is easy to check directly that the generating 
function obeys the following recursion relation 



G t +i = («aAb + ?^bAc + "-cAd) 
+ n D (pi A D G t + P2^cGt + P3AbG( 



(16) 
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In the limit of large avalanches we believe that this 
generating function tends to some universal function 
Goo(Aa, Ab, Ac, Ad)- Thus, ( |l6|) becomes a closed equa- 
tion for this function. This equation depends on two set 
of parameters n and p. In the critical state, however, 
these are not independent. The relation between them 
can easily be found if we note that "below" the critical 
point, when all avalanches are finite, the average num- 
ber of species A, B, C and D in the toppled sites can be 
simply related to the generating function 
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(17b) 
(17c) 
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As expected, these numbers become divergent at the crit- 
ical state. Calculating the concentrations of species A, 
B, C and D after the avalanche and equating them to the 
concentrations before the avalanche 



(Na) 



(Na) + (N B ) + {Nc) + (N D ) 



= n A , ■ 



(18) 



we will reproduce again the relations (|^), thus proving 
the self-consistency of this mean-field description. 

It is hardly possible to solve the forth order equation 
(|l|) explicitly. We can simplify it further (still within 
the mean-field approximation) if we take into account 
the fact that the average number of particles leaving the 
site on toppling is equal to p and each generates on av- 
erage equal sub-avalanches. Then, our equation for the 
universal generating function becomes 



where we put Aa = Ab = Ac = Ad = A. A similar re- 
cursion relation naturally appears in any mean-field de- 
scription of the avalanche process ||-^| . 

In the simplest case p = 2 the solution of equation ( |l9| ) 
can be written in the closed form 



Goo (A) = 



1- v /l-4A 2 np(l-np) 
2Atid 



(20) 



which, after expanding it as a series in A, leads to the well 
known mean-field probability distribution for the sizes of 
avalanches at the critical state 



P(N) ~ N~ 3 / 2 . 



(21) 



Actually, this asymptotic behavior does not depend on 
the value of p and the mean- field exponent 3/2 is univer- 
sal for any critical branching process pT| . 



III. SPANNING TREES IN THE SANDPILE 
MODEL 

A. Basic Properties 

Now let us go back to the original formulation of the 
sandpile model. It has been shown by Dhar |l^] that this 
model is actually exactly solvable. 

To derive these results it is convenient to reformulate 
the dynamical rules of the sandpile model as follows. We 
consider the model on a square lattice C of N sites la- 
beled by 1, 2, N . Each boundary site is connected by 
a bond to the additional site 7k- (the root) which plays 
the role of the sink for the particles. This auxiliary site, 
although unphysical, will be very convenient for all our 
further constructions. The N x N matrix of the discrete 
Laplacian Ay has non-zero diagonal elements equal 
to the number of neighboring sites of i and non-zero off- 
diagonal elements Ay = — 1 for all pairs of adjacent sites 
i and j. The addition of sand corresponds to increasing 
the height of the pile by unity at a site of the lattice 
chosen at random (except *). If the height at any site i 
exceeds its critical value A# , that site topples and all the 
variables Zj, (j = 1, N) are updated according to the 
rule 



Goo(A) = A((l - n D ) + n D [Goo(A)f ), 



(19) 



The process stops when all the heights of the resulting 
configuration ip — {zi} do not exceed their critical values. 
Such a configuration is called a stable configuration. 

Yet, it is not clear whether the dynamics of the model 
is well defined because during the toppling process a con- 
figuration may occur with two or more unstable sites. We 
have to make sure that the resulting stable configuration 
does not depend on the order of their topplings. This can 
easily be verified if we note that after the toppling of two 
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unstable sites i and j in an arbitrary order one gets a con- 
figuration in which height z% at any site k of the lattice 
C decreases by + Aj/.. This expression is invariant 
under the exchange of i and j. By a repeated use of this 
argument, one obtains that after any avalanche the same 
final stable configuration is produced irrespectively of the 
order of topplings of unstable sites |l2| . 

Another possible ambiguity is due to the fact that the 
updating procedure may in principle enter a nontrivial in- 
finite cycle. Nevertheless, one can easily prove that this is 
impossible by noting that topplings in the interior of the 
lattice does not change the total amount of sand in the 
system and every toppling on the boundary decreases this 
value. Here, no cycle can have topplings at the bound- 
ary. Next, the sand on the boundary will monotonically 
increase if there is any toppling one site away. This can- 
not happen in the infinite cycle, thus there can be no 
topplings one site away from the edges. By induction, 
there can be no topplings at an arbitrary distance from 
the boundary, thus, there can be no infinite cycle 

The Markovian evolution of the model implies that the 
set of all stable configurations S can be divided into two 
subsets: recurrent R and transient S \ R. By definition, 
the first subset includes those stable configurations that 
can be reached from all others by sequential addition of 
particles. This subset is not empty because there is at 
least one such a configuration -00 — { z i — A^} with all 
the heights equal to their critical values. Then, the very 
general result of the theory of Markovian chains states 
that once the system gets into the set of reccurent con- 
figurations, it will never get out under the further evolu- 
tion. All non-recurrent stable configurations are usually 
called transient. 

To calculate the number of recurrent configurations 
jHI , let us consider the space of all possible (including 
unstable) configurations obtainable from fa by addition 
of particles. Then, we define two such configurations {zi} 
and {z[} as equivalent if and only if under topplings they 
evolve to the same stable configuration. This means that 
there exists some integers rj such that 

N 

z[ = Zj — rjAjj, for all i. (23) 

Thus, if we represent configurations {z{\ by points in a 
TV-dimensional hypercubical lattice with basis vectors ej, 
the set of equivalent points forms a super-lattice with 
basis vectors 

JV 

Ei=^2 A u" g j"' i = 1 to JV. (24) 

Since every class of equivalent configurations corresponds 
to some unique recurrent configuration, the volume of the 
unit cell of the super-lattice must be equal to the number 
of distinct recurrent configurations. Thus one gets 



Af R = detA. (25) 

It is important to note that, although the geometri- 
cal shape of the set R in this TV-dimensional space is 
quite nontrivial, the copies of R can be arranged to give 
a simple periodic tilling of the space. In other words, 
the set R of all reccurent configurations has the topology 
of a A-dimensional torus. Therefore, the addition of a 
particle to an arbitrary site i of the lattice can be rep- 
resented by a jump of the point on this torus along the 
corresponding unit vector (in the positive direction). 
Then, if one drops particles onto different sites of the 
lattice with equal probabilities, one forces the represent- 
ing point on the torus to move randomly (but always in 
the positive directions of the axes) from one site of the 
torus to another. Once this random motion covers the 
torus uniformly, all the reccurent configurations will ap- 
pear in the Markovian evolution of the system with equal 
probability. 

B. Forbidden Sub-Configurations and Spanning 
Trees 

To describe in detail the set of recurrent configurations, 
we introduce first the important concept of forbidden sub- 
configuration |l4f| . It is defined as an arbitrary subset of 
sites of the lattice T C C such that all its corresponding 
height variables {zj}, j E J-, satisfy the inequalities 

Zi < deg^(i) = (~ A y)> foralHE-F, (26) 

where we denote by degyr(i) the number of bonds con- 
necting the site i with the other sites of the subset T . 

Any stable height configuration that contains no for- 
bidden sub-configurations is called an allowed configu- 
ration. We will prove soon that the set of all allowed 
configurations A is the same as the set of all reccurent 
configurations R. 

First, let us show that the set of allowed configurations 
is closed under the dynamics of the sandpile model |L2l . 
In other words, we want to show that once the system 
gets into the set of allowed configurations it will never get 
out. Assume the contrary. We then note that the addi- 
tion of particles only increases heights and, hence, cannot 
create a forbidden sub-configuration. This means that 
there exists an allowed configuration ip such that after 
a single toppling it becomes the configuration ip 1 which 
contains a forbidden sub-configuration T . Suppose the 
toppling occurs at site % E T ' . From the toppling rules 
and the definition of forbidden sub-configuration one gets 
that if T is a forbidden sub-configuration in ip', then the 
set !F\i (obtained from T by deleting the site i) is forbid- 
den in the initial configuration ip. This contradicts our 
assumption that ip is allowed. Finally, since the recur- 
rent configuration fa is allowed and all other recurrent 
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configurations are reachable from this particular one, it 
follows that all recurrent configurations are allowed 

RC A and Af R < Ma- (27) 

Another important notion is that of toppling from the 
root. It gives a simple recursive procedure determining 
whether a given stable configuration -0 is allowed. 

To organize such a toppling, we topple the auxiliary 
site * as if it were an ordinary site of the lattice. We 
could equivalently define it in terms of the simultaneous 
dropping of particles onto each boundary site of the lat- 
tice (i.e. one particle at each site on the edges and two 
particles onto each corner site). If no boundary sites be- 
come unstable then the original configuration tp is not 
allowed. Using the definition one can easily show that in 
this case the lattice itself plays the role of forbidden sub- 
configuration: T — C. More generally, if the avalanche 
initiated by the toppling from the boundary has stopped 
after a few topplings and some set T of lattice sites re- 
mains untoppled, then this set T is nothing but the for- 
bidden sub-configuration associated to the initial config- 
uration ip. 

In such a toppling no lattice sites can topple more than 
once |Q. To prove this, let us assume that some site i 6 
C has toppled for the first time after its nearest neighbor 
j. Then, i would topple for the second time only after 
the topplings at all its neighbors, including j. Therefore, 
to topple i twice, we have first to topple the neighboring 
site j twice. Using the same arguments for j and for sites 
that toppled earlier than j, we conclude that to topple i 
more than once, we have to topple some boundary sites 
a second time. This is obviously impossible due to the 
loss of particles after the first toppling of boundary sites. 

Thus, we have proved that the initial configuration ip is 
allowed if and only if the toppling from the root generates 
an avalanche process under which each site of the lattice 
topples only once and, hence, the height configuration 
remains unchanged. If the initial configuration is not 
allowed, then some untoppled sites will survive and will 
form a forbidden sub-configuration. 

It is possible to visualize the avalanche process initi- 
ated at the root. To this end let us consider the dynam- 
ical process step by step. Initially only some of the sites 
at the boundary can become unstable. They will topple 
at the moment t = 1 and transfer particles to neigh- 
bors that can become unstable in their turn. Similarly, 
all sites unstable at the moment t topple simultaneously 
and produce new unstable sites that will topple at the 
next time step t + 1. Consider an arbitrary site i. Let 
t be the time step at which this site becomes unstable 
and t + 1 the moment at which it topples. This means 
that at least one of its nearest neighbors toppled at the 
time step t. Let £ be the number of such neighbors that 
toppled at the moment t. The height z, should obey the 
following inequalities 



A« - f < Zi < A u , (28) 

because otherwise it could not topple at the moment t+1 
as was assumed. Now, if £ = 1, we simply mark in red the 
only bond connecting the site i with the neighboring site 
that toppled at the earlier moment t. In the other case, 
when £ > 1, we select from these £ possibilities only one 
bond to mark in red, dependent on the height Zi at site i. 
To this end we create an order of preferences by enumer- 
ating all bonds incident to the site i in an arbitrary but 
fixed manner. We then choose from the £ candidates the 
bond that occupies the [z% + £ — A^)-th position in this 
list of preferences. For example; if the northward bond of 
the site i has been allocated the number 1; the eastward 
2; the southward 3; and the westward 4 and if the two 
neighboring sites of i at the north and at the south topple 
at the time step t, then if Zi = 4, we mark the southward 
bond (numbered by 3) red as it is the greater of the two 
in the list of preferences. This algorithm makes it possi- 
ble to avoid any ambiguity in the choice of red bonds and 
the construction of a unique graphical representation of 
a given allowed configuration. 

As all sites of the allowed configuration must topple, 
the graph 7^ formed by the red bonds will cover all sites 
of the lattice. Such a graph is usually called the spanning 
graph. The fact that each site of the lattice has toppled 
only once implies that the spanning graph T+ contains no 
loops. Such a graph is called a spanning tree. A spanning 
tree having one site (the root ★) distinguished from all 
others is called a rooted spanning tree. Since there is 
only one path from any site of the rooted spanning tree 
to the root, we may uniquely orient this path so that 
all of its red bonds will be supplied with arrows in the 
direction of the root |l7j . 

It is not difficult to check that, given a rooted spanning 
tree on the lattice, one can easily reconstruct the height 
configuration using the same list of preferences as above. 

Thus, starting with the definition of forbidden sub- 
configurations, we have proved the one-to-one correspon- 
dence between allowed sandpile configurations on the lat- 
tice and rooted spanning trees on the same lattice @. 

C. Correlation Functions and the Continuous Limit 

The spanning tree is actually a strongly correlated 
object. This explains why the absolutely uncorrelated 
heights of an arbitrary initial configuration become cor- 
related in the recurrent state during the course of relax- 
ation. The tree analogy provides a useful representation 
for the determination of the statistical properties of the 
sandpile model in the critical state. The tool used to 
enumerate rooted spanning trees is given by the famous 

Kirchhoff 's theorem jl?]] . If to any bond of the lattice 
C, whose adjacent sites i and j are different from the 
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root, the weight X ij IS ascribed, then the determinant of 
the matrix 



G(0,1) -G(0,0) 



'Y,k x ik , if i = 3 
&>ij{x) = { -Xij , if |* — j\ = 1 
, otherwise. 



is a partition function of the rooted spanning trees on the 
lattice. 

The fact that we treat here the weights i.y and Xji as 
different variables implies that we consider the arrow out- 
going from site i and directed towards j as different from 
its opposite. When — 1 for all adjacent sites i and j 
this matrix coincides with the discrete Laplacian and its 
determinant gives the total number of rooted spanning 
trees. Hence the number of distinct allowed configura- 
tions is given by 



N A = det A. 



(30) 



The equivalence between the set A of all allowed configu- 
rations and the set R of all reccurent configurations now 
follows immediately from Eqs.( p5| , ^7||3C| ). 

Using the KirchhofPs theorem one can determine, 
in principle, all of the correlations between different 
branches of spanning trees. Although the method of such 
a calculation is far from being novel, it is still worth re- 
calling the principal ideas. 

Any modification of the weights of a finite number of 
lattice bonds is called a local defect of the lattice. For 
example deleting the bonds or inserting additional ones 
can be considered as a proper local defect. The difference 
between a discrete Laplacian of the new lattice A' and 
that of the initial one A is referred to as the defect matrix 
(5. The locality condition simply implies that only a finite 
number of the rows and columns of the defect matrix S 
have non-zero elements. It follows that the the ratio of 
the number of spanning trees on the new lattice £ to the 
number of all trees on the original one C is given by an 
easily calculable determinant 



Prob((5) 



det A' 
det A 



= det(l + G<5) 



(31) 



where 1 is the unit matrix and G = A -1 is the lat- 
tice Green's function. On the planar square lattice the 
Green's function only depends on the distance between 
the sites i and j and has the following integral represen- 
tation 

If cosna cosm/3-1 

G 7i, to) - G O, 0) = — / — — dadp. 

87H J 2 — cos a — cos p 


(32) 

Here n and to are the numbers of columns and rows be- 
tween the sites i and j respectively. This integral can be 
calculated explicitly with the results 



G(1,0)-G(0,0) = --, (33) 



(29) G (»-»)- G ( -»)-K 1 + i + 5 + - + 2^ L r) 



(34) 

Using these values and the discrete Laplace equation as a 
recursion relation, one can find all other elements of the 
matrix Gy. 

For example, let us find the probability of having the 
height Zi = 1 at some site i deep within the lattice jl9). 
One can notice that decreasing the height at the site i 
by 1 — so that its height becomes equal to — one 
ends up the forbidden sub-configuration that consists of 
only the site i itself. In the language of spanning trees 
this corresponds to those trees that cover site i by a leaf 
bond (by deleting the bond site i becomes disconnected 
from the rest of the tree). Due to the equivalence of the 
4 positions of the leaf bond, the unit height probability 
we are interested in can be expressed as 



Prob(zi 



1) 



1 Af' _ 1 det A' 
4 If ~ 4 det A ' 



(35) 



Here A/"' is the number of spanning trees on the new lat- 
tice £ where all the bonds connecting the site i with its 
four nearest neighbors are deleted and, instead, the site i 
is directly connected to the root *. Labeling the nearest 
neighbors of the site i in a clockwise fashion as N, E, S 
and W, we can write the corresponding defect matrix 



(36) 



Direct evaluation of the determinant Eq.(j3l|) then leads 
to the result 







i 


N 


E 


s 


w 


i 


( 


-3 


1 


1 


1 


1 \ 


N 




1 


-f 











5= E 




1 





-1 








S 




f 








-1 





W 


V 


f 











-v 



2 / 2 
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0.073 636. 



(37) 



The calculation of all other height probabilities requires 
more sophisticated ideas and can be found in pcfl . 

Similarly, we could calculate the asymptotic probabil- 
ity that sites i and j (both deep within the lattice) sep- 
arated by the distance r will both have height 1 in the 
reccurent state. It is more instructive however to em- 
ploy another approach to find the asymptotics of this 
two-point correlation function. 

Namely, we can reinterpret the partition function of 
the spanning trees (given by KirchhofPs theorem) as be- 
ing the partition function of some artificial statistical sys- 
tem. To this end we define at each site i of the lattice £ 
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the pair of anticommuting variables 0, and 8* (its con- 
jugate). Then, using Berezin's definition of the integral 
over anticommuting variable s pjj we can rewrite the de- 
terminant of the matrix Eq. (|29[) as 



Z = detA i3 = J d6>i...<20^exp^0*A%. 



(38) 



In the continuous limit this partition function defines a 
conformal field theory with the central charge c = — 2 
and the Hamiltonian (after integration by parts) 



H 



3^8*3^9 d 2 r. 



(39) 



The Green function of the field 9 and its conjugate 9* co- 
incides with the asymptotic behavior of the lattice Green 
functions Ga 



! (r x )0(r 2 )) = G(n - r 2 ) = G(0) - ±- In |n - r 2 | 



(40) 



Changing the weight of a given directed bond in Eq.(38), 
say Xij , and then again taking the continuous limit we 
will get the current operator corresponding to the fixed 
arrow on the spanning tree 



= 9*d u e. 



(41) 



Hence we can immediately determine the asymptotic be- 
havior of the arrow-arrow correlation function on the 
spanning tree 



0»^(o)> = 0m)(>) 



1 X nXy 



At: 2 



(42) 



Similarly, using the defect matrix 5 
the local energy operator 



we can define 



8^9*8^9 



and calculate its two-point correlation function 
( £( r) £ (0)> = ( £ ) 2 - i ^A 



(43) 



(44) 



This formula gives the asymptotic behavior of the two- 
point correlation function of unit heights in the sandpile 
model |H| (up to proper normalization of the local en- 
ergy operator in the continuous limit). Moreover, on the 
boundary of the lattice only this operator determines the 
correlations of all other heights p3,p3| . 



D. Waves of topplings 

The study of avalanches requires a further extension 
of the spanning tree representation. In order to make it 



possible we have to reorganize the topplings inside the 
avalanche into successive waves of topplings |l5j . As has 
already been mentioned, the dynamics of the sandpile 
model admits an arbitrary order of topplings of unstable 
sites during an avalanche. We choose a particular — but 
quite natural — order amongst all others. Namely, let 
us drop a particle onto the critical site i in an allowed 
configuration ip. We topple it once and then topple all 
sites that become unstable keeping the initial site i out 
of the second toppling. We call the set of sites toppled in 
this way the first wave of topplings, denote it T\ . After 
the first wave has gone out we topple the site i a second 
time and continue the avalanche, never permitting this 
site to topple a third time. The set of relaxed sites in the 
period after the first wave is called the second wave of 
topplings, denote it J-2- The process continues produc- 
ing intermediate configurations ipi,ip2> until the site 
i undergoes the maximum number of topplings and the 
avalanche stops. 

In complete analogy with the case of toppling from the 
root, the sites of any set Tk topple only once during the 
£;-th wave. The set Tk has no holes: otherwise the subset 
of sites corresponding to the hole would form a forbidden 
sub-configuration of the initial recurrent configuration ip. 
The compact cluster Tk of sites toppled in the &-th wave 
forms the forbidden sub-configuration Tk of the configu- 
ration ipk- Indeed, the inequality (26) holds for the height 



of every site in the cluster provided that all other sites of 
the lattice remain untoppled during that wave. 

Thus, the avalanche starting at i is represented as a 
collection of T waves, T > 1, which we also denote by 
Fk, (1 < k < T). All waves with numbers 1, ...,T - 1 
have the site i strictly inside Tk- The cluster Tt of the 
last wave, T, has the site i on its boundary. Indeed, the 
avalanche can stop only if the site i has at least one neigh- 
boring site not belonging to the last wave and therefore 
giving no contribution to the growth of the height Zi in 
ipr- 

To find the tree representation of waves JTj| , we con- 
sider the sandpile model on a new lattice that consists 
of the original one C and the additional bond (*i) con- 
necting a given site i with the root Accordingly, we 
have to change one element of the matrix of the discrete 
Laplacian which determine the toppling rules 



A, 



A 4l + 1, 



(45) 



Now, to construct the spanning tree corresponding to 
the wave of topplings, let us consider again the toppling 
from the root *. In this toppling, particles drop onto the 
boundary sites of the lattice and one particle goes along 
the bond connecting the root * with the origin of the 
avalanche i. Let us consider these topplings separately 
(again choosing a particular order of the topplings). We 
start by sending one particle to the site i. The avalanche 
starting at the site i spreads over some part T\ of the 
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lattice. The site i can topple only once (as it is always 
connected to the root in this new lattice) hence T\ is the 
first wave of topplings. We drop then all other particles 
onto the boundary sites causing another avalanche that 
should cover the rest of the lattice C\T\. This avalanche 
returns the configuration updated after the first wave to 
the initial state and can be termed the inverse wave. 

We can generalize the procedure described above. Let 
us perform k topplings from the root on the new lattice 
with the additional bond and allow k particles to 
pass through the bond They will initiate exactly 

k waves T\, having i as origin. If we then drop 
particles once onto the boundary sites, thus causing an 
inverse wave, we will return the configuration of heights 
to the previous one (V>fc-i), corresponding to the (k — 1)- 
th wave. This means that there is a complete commuta- 
tivity between waves and inverse waves and they can be 
performed in an arbitrary order. Moreover, starting from 
a recurrent configuration on the new lattice, by deleting 
one particle at site i and initiating only inverse waves, 
one can reverse completely the evolution of the sandpilc 
model. 

We have seen in the previous sections that the top- 
pling process naturally draws spanning trees on the lat- 
tice. Applying the algorithm described above to the re- 
current configurations on our new lattice with the addi- 
tional bond (★£), we get a new set of spanning trees. All 
spanning trees now can be divided into two sets. The 
first one consists of only those spanning trees that do 
not contain the bond and, therefore, on deleting the 
bond coincides with the set of one-rooted spanning trees 
%, on the old lattice C. The second set consists of those 
trees that contain the bond On deleting this bond 
(thus returning to the original lattice C) a subtree % gets 
disconnected. Considering the site i as the second root of 
component % we obtain a two-rooted situation, where a 
spanning tree on the lattice consists of two disconnected 
clusters % and T+. 

Thus, in addition to the one-to-one correspondence be- 
tween the recurrent states and one-rooted spanning trees, 
we get the one-to-one correspondence between all waves 
of topplings and all two-rooted spanning trees. The 
avalanche is displayed now as a collection of successive 
two-rooted trees. 



E. Waves and Green functions 

The graph representation of waves enable us to link 
the toppling process and the lattice Green function Gij . 

Theorem [ [l5| . For an arbitrary lattice C with the root 
*, the Green's function is given by the ratio 



— 



(46) 



where is the number of two-rooted spanning trees 

having the roots * and i, such that both vertices i and j 
belong to the same one-rooted cluster and Af is the total 
number of one-rooted spanning trees on C. 

As a simple consequence of the theorem we conclude 
that Gij is equal to the expected number of topplings 
at the site j during the avalanche caused by adding a 
particle at site i fH| . Indeed, as each wave corresponds to 
only one toppling of its sites, then the expected number 
of topplings should coincide with the expected number 
of waves involving the site j and, hence, should be equal 
to Gij. 

Yet another result is necessary to derive the proba- 
bility distributions of waves. Namely, as is known from 
conformal field theory arguments |24-26|, the area s and 
the perimeter I of a finite cluster of two-rooted spanning 
tree (or any wave of topplings) are related to its radius r 
as 



r 2 and I 



r 5/4 



(47) 



To find the size distribution of general waves |l5| , we 
consider all possible waves that belong to the avalanches 
starting at a fixed point i deep within the lattice. All 
these waves are in in one-to-one correspondence with the 
recurrent configurations of the sandpilc model on the lat- 
tice with an additional bond connecting the sites * and i. 
As all recurrent configurations have the same probability 
it follows that all general waves are also equally likely. 

Let us now consider waves on the lattice initiated at 
site i without any reference to the particular avalanches 
they belong to. In other words, we consider a sequence of 
avalanches initiated at the same site i. These avalanches 
consist of waves. We ignore the pauses between the 
avalanches and study the properties of the collection of 
waves in general. Using the above theorem we can esti- 
mate the probability that the radius r of the wave is not 
less than the distance between the sites i and j as 



Prob(r > \i — j\) ~ Gij ~ In \i — j\ 
The corresponding probability density 

dr 

P w {r) dr ~ — , 



(48) 



(49) 



is only normalizable provided that the size of the lattice L 
(upper cutoff) and its spacing a (lower cutoff) are fixed. 
Using scaling relations Eq.(f47|) and 



P w (s) ds~P w (l) dl~P w {r) dr, 



(50) 



we can rewrite this probability distribution either in 
terms of the area s or the perimeter I of the waves of 
topplings. 

We have seen above that the last wave corresponds to 
the rooted subtree having its root at the boundary. The 
alternative way to get its tree representation is to cut a 
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bond and disconnect a branch from the one-rooted tree 
representing the recurrent configuration that appears af- 
ter the last wave. 

A general wave differs from the last one only in the 
location of the root. The position of the root in the gen- 
eral wave is distributed uniformly over the whole area, in 
contrast to its location at the boundary for the last wave. 
Therefore, given a general wave of area s and perimeter 
I, an additional factor l/s ~ s -3 / 8 appears in the prob- 
ability distribution of the last waves. Hence, we get the 
size probability distribution of the last waves p7[ 

Pi(s) ds ~ P w (s) s- 3/8 ds ~ s~ 11/8 ds. (51) 

The critical exponents of the general waves are not 
related directly to the exponents of the avalanche dis- 
tributions because by considering the general wave we 
lose information as to the concrete avalanche to which 
it belongs. Fortunately, there exists a situation where 
all avalanches consist of one and only one wave and the 
avalanche distribution coincides with the distribution of 
general waves. Namely, avalanches starting at the bound- 
ary consist of only one wave because the second toppling 
is impossible due to loss of particles. 

Let us consider the case of an infinite wedge, with in- 
ternal angle 7, related to the half plane by the conformal 
mapping w — z 7 ^. It is known from the theory of com- 
plex variables that the Green function of the Laplacian 
in the region inside the wedge has the form 



G{z,z) ~ 1m(z-* / ' y ) 



(52) 



where z and z,are the complex coordinates on the plane. 
Thus, the function G(r) decays as 



G(r) 



-7T/7 



(53) 



for all directions save the arms of the angle. This leads 
to the distribution 



P b (r) dr ~ r- 1 -"^ dr 



(54) 



Again, using the relations s ~ r 2 and Pb(s) ds ~ Pb(r) dr, 
we can rewrite this distribution in the form 



Pb(s) ds 



ds 



(55) 



which corresponds to the critical exponent r = 1 +71-/27 
for the area distribution of boundary avalanches (2f| . 

The angle 2tt is of special interest. In this case 
avalanches start at the top of a cut on the plane. The 
geometry of the avalanches closely resembles the one oc- 
curring deep within the lattice. So one could expect that 
the corresponding critical exponent r = 5/4 should be 
close to that occurring in the bulk of the lattice. 



IV. AVALANCHE DISTRIBUTION EXPONENTS 

Let R denote the linear extent (diameter of the 
avalanche); S, the number of distinct sites where at least 
one toppling occurs; V, the total number of topplings in 
the avalanche. We denote the number of waves of top- 
plings it consists of as T. All these quantities are random 
variables. It is usually assumed that their probability 
distributions have the following power-law asymptotics 
in the thermodynamic limit 



P a (R) dR ~ iT TR dR 
P a {S) dS ~ S~ TS dS 
P a (V) dV ~ V~ TV dV 
P a {T) dT ~ T~ TT dT 



(56) 
(57) 
(58) 
(59) 



The exponents tr,ts,tv and tt can be related to each 
other if we make the simple scaling assumptions that 
R, S, V and T scale as some powers of each other, i.e. 
if we assume that for all avalanches with a fixed value 
of one variable (say R) the other variables (S, V and 
T) have strongly peaked conditional probability distri- 
butions Q . 

Since every avalanche cluster is nothing but a superpo- 
sition of compact clusters of waves of topplings, we may 
speculate that the avalanche clusters are also compact 



S ~R 2 



(60) 



The total number of topplings in the avalanche can then 
be related to its diameter via some new exponent y 



V ~ R z 



where y > 0. 



(61) 



Finally, from the fact that in each wave of topplings the 
sites of the lattice topple only once, we conclude that the 
same exponent should appear in the relation between the 
number of waves in the avalanche and its linear extent 



T ~ R v . 



(62) 



The probability distribution P a {T) of the number of 
waves in the avalanche should be consistent with these 
scaling relations. It should also lead to the logarithmic 
dependence of the average number of waves on the lat- 
tice size, (T) ~ InL. The only distribution that meets 
all these assumptions is 



P a (T)dT~^. 



(63) 



Now, using scaling relations Eq. ( |60y61|j62| ) we can express 
the critical exponents of all other probability distribu- 
tions in terms of the only unknown scaling exponent y 



TR = l + y, 

t s = 1 + y/2, 

T V = (2 + 2y)/{2 + y), 

T T = 2. 



(64) 
(65) 
(66) 
(67) 
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In the particular case y = these exponents coincide 
with those of waves of topplings. In the opposite case, 
y = 1, they reproduce the mean-field critical exponents. 
Thus, we may expect that the probability distribution 
of avalanches should corresponds to some intermediate 
value < y < 1. 

The more detailed description of the avalanche pro- 
cess needs new assumptions about the average space- 
time structure of the avalanche process and its decom- 
position into the waves of topplings. On the basis of 
numerical simulations it can be assumed that, typically, 
the avalanche process consists of two regimes. In the 
first regime one witnesses a relatively fast growth of 
the avalanche with each subsequent wave escaping the 
boundary of its predecessor, while in the second regime 
one observes a relatively slow squeezing of waves towards 
the last wave of topplings that touches the initial site. 
Assuming that in the thermodynamic limit the second 
stage completely dominates the first one can expect that 
the area decrement Ast between two successive waves St 
and St+i scales along with their average size Sk as 



As t ~ a?. 



(68) 



This scaling law is consistent with all our earlier assump- 
tions. The exponent a can be related to the scaling ex- 
ponent y that determine the avalanche probability distri- 
butions. Indeed, we have 



As t ~ s?At, and Ar t 



r? a_1 Ai, 



(69) 



where t is the number of wave (At — 1 for two successive 
waves). If we integrate these relations up to the size of 
the avalanche R and the total number of waves T, we 
reproduce Eq.(|62|) with the scaling exponent y — 2 — 2a. 
The attempt to find a from more complicated scaling 
arguments Q gave a = 3/4 and, hence, y = 1/2. 

Although these critical exponents are in agreement 
with the results of the real space renormalization group 
calculations ^,||], it is worth mentioning again that the 
above construction is based only on some self-consistent 
hypothesis. It is not clear at the moment how they are 
supported by the numerical simulations. 

An alternative hypothesis about inner structure of the 
avalanche process assume the existence of the universal 
asymptotics of the probability distribution of subsequent 
waves for a given size of the preceding wave, P(st+i\st) 
pl[ . Namely, on the basis of numerical simulations it 
was suggested that 



P(st+i\s t ) 



't+1 
i 7 

't+1 



if s t+ i < s t 
if s t +i > s t 



(70) 



with (3 = 3/4 and 7 = 5/4. The derivation of this scal- 
ing laws, if they exist, requires further analysis of the 
avalanche process. 



We would like to thank Deepak Dhar for many dis- 
cussions and correspondence and James Harris for useful 
comments on the manuscript. EVI is grateful to Prof. 
Chin-Kun Hu for kind invitation to make this talk at the 
conference "StatPhys-Taipei-1997" . VBP thanks Dublin 
Institute for Advanced Studies for kind hospitality. 



[1] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59 
(1987) 381. 

[2] A.N. Kolmogorov, C. R. Acad. Sci. USSR 30 (1941) 301. 
[3] M. Kardar, G. Parisi and Y.-C. Zhang, Phys. Rev. Lett. 
56,889 (1986). 

[4] LP. Kadanoff, S.R. Nagel, L. Wu and S.M. Zhou, Phys. 

Rev. A 39 (1989) 6524. 
[5] E.V. Ivashkevich, Phys. Rev. Lett. 76 (1996) 3368. 
[6] C. Tang and P. Bak, Phys. Rev. Lett. 60 (1988) 2347. 
[7] D. Dhar and S.N. Majumdar, J. Phys. A 23 (1990) 433. 
[8] S.A. Janowsky and C.A. Laberge, J. Phys. A 26 (1993) 

L973. 

[9] K. Christensen and Z. Olami, Phys. Rev. E 48 (1993) 
3361. 

[10] Y.-C. Zhang, Phys. Rev. Lett. 63 (1989) 470. 
[11] M.E. Fisher and J.W. Essam, J. Math. Phys. 2 (1961) 
609. 

[12] D. Dhar, Phys. Rev. Lett. 64 (1990) 1613. 

[13] M. Creutz, Comput. Phys. 5 (1991) 198. 

[14] S.N. Majumdar and D. Dhar, Physica A 185 (1992) 129. 

[15] E.V. Ivashkevich, D.V. Ktitarev and V.B. Priezzhev, 

Physica A 209 (1994) 347. 
[16] CM. Fortuin and PW. Kasteleyn, Physics 57 (1972) 536. 
[17] F. Harary, Graph Theory (Reading, MA: Addison Wesley, 

1990). 

[18] F. Spitzer, Principles of Random Walk (Van Nostrand, 
New York, 1964). 

[19] S.N. Majumdar and D. Dhar, J. Phys. A 24 (1991) L357. 

[20] V.B. Priezzhev, J. Stat. Phys. 74 (1994) 955. 

[21] J. Zinn- Justin, Field Theory and Critical Phenomena 
(Oxford, Clarendon, 1993). 

[22] E.V. Ivashkevich, J. Phys. A 27 (1994) 3643. 

[23] J.G. Brankov, E.V. Ivashkevich and V.B. Priezzhev, 
J. Phys. I France 3 (1993) 1729. 

[24] J.L. Cardy,P/iase Transitions and Critical Phenomena 
vol 11, eds C Domb and J L Lebowitz (London: Aca- 
demic Press, 1987) 

[25] H. Saleur and B. Duplantier, Phys. Rev. Lett. 58 (1987) 
2325. 

[26] A. Coniglio, Phys. Rev. Lett. 62 (1989) 3054. 

[27] D. Dhar and S.S. Manna, Phys. Rev. E 49 (1994) 2684. 

[28] E.V. Ivashkevich, D.V. Ktitarev and V.B. Priezzhev, J. 

Phys. A 27 (1994) L585. 
[29] L. Pietronero, A. Vespignani and S. Zapperi, Phys. Rev. 

Lett. 72 (1994) 1690. 
[30] V.B. Priezzhev, D.V. Ktitarev and E.V. Ivashkevich, 

Phys. Rev. Lett. 76 (1996) 2093. 
[31] M. Paszuski and S. Boettcher, Phys. Rev. E (to appear). 



10 



